import numpy as np
from scipy.signal import czt, find_peaks, filtfilt
# 控制干扰峰个数的 幅值判断
interference = 0.001

# 控制通用波峰检测的 幅值判断
general_peak_threshold = 0.1


def determine_type(variable):
    """类型判断工具"""
    if isinstance(variable, list):
        return "list"
    elif isinstance(variable, np.ndarray):
        return "numpy.ndarray"
    elif isinstance(variable, (int, float)):
        return "number"
    elif isinstance(variable, str):
        return "string"
    elif isinstance(variable, dict):
        return "dict"
    elif callable(variable):
        return "callable"
    else:
        return "unknown"
#信号预处理模块
class SignalPreprocessor:
    """信号预处理模块"""
    interference = 0.001  # 控制干扰峰判断的幅值阈值

    HF_a = [1, -3.590148110152597, 4.873443614170354, -2.963870062130241, 0.681973883426416]
    HF_b = [0.775707823645609, -3.094022833216837, 4.636644012395594, -3.094022833216837, 0.775707823645609]
    LF_a = [1, -2.174805506915435, 2.357402921272809, -1.282544135294374, 0.306064117903952]
    LF_b = [0.036936953516246, 0.031977931425171, 0.061290122632184, 0.031977931425171, 0.036936953516246]
    PF_b = [0.142537860520234, -0.215022543769845, 0.149319417312823, -0.215022543769845, 0.142537860520234]
    PF_a = [1, -2.67035772773509, 2.94351262911759, -1.63009314766181, 0.400438754415306]

    @staticmethod
    def ppg_process(channle_data, fs):
        """
            计算PPG,对信号进行处理
        """
        # 上面那个是对PPG信号做异常处理的，下面这个是做平滑的。
        abnormal_remove_ppg_channle = SignalPreprocessor.abnormal_remove(channle_data)
        abnormal_remove_ppg_channle = SignalPreprocessor.smooth_filt(abnormal_remove_ppg_channle, 6)

        filt_process_ppg_channle = SignalPreprocessor.filt_process(abnormal_remove_ppg_channle, fs)

        # 这个做动作干扰滤波的：
        motion_filt_process_ppg_channle = SignalPreprocessor.motion_filt_process(filt_process_ppg_channle)

        return motion_filt_process_ppg_channle

    @staticmethod
    def acc_process(channle_data, fs):
        """
            计算ACC,对信号进行处理
        """
        # 上面那个是对ACC信号做异常处理的，下面这个是做平滑的。
        abnormal_remove_acc_channle = SignalPreprocessor.smooth_filt(channle_data, 6)

        filt_process_acc_channle = SignalPreprocessor.filt_process(abnormal_remove_acc_channle, fs)
        # 这个做动作干扰滤波的：
        motion_filt_process_acc_channle = SignalPreprocessor.motion_filt_process(filt_process_acc_channle)
        return motion_filt_process_acc_channle

    @staticmethod
    def abnormal_remove(channle_data):
        """去除信号异常值"""
        """
            去除异常值
            """
        HF_a = [1, -3.590148110152597, 4.873443614170354, -2.963870062130241, 0.681973883426416]
        HF_b = [0.775707823645609, -3.094022833216837, 4.636644012395594, -3.094022833216837, 0.775707823645609]
        # 去均值
        channle_data = channle_data - np.mean(channle_data)
        # 一阶差分
        diff_channel_data = np.diff(channle_data)
        diff_channel_data = np.insert(diff_channel_data, 0, 0)
        # 找一阶差分波谷
        valley_index_list, _ = find_peaks(-diff_channel_data, distance=6)
        # 相邻两个波谷之间最大值位置找到被定义为波峰位置
        peak_index_list = []
        for temp_index in range(len(valley_index_list) - 1):
            peak_index = np.argmax(diff_channel_data[valley_index_list[temp_index]: valley_index_list[temp_index + 1]])
            peak_index += valley_index_list[temp_index]
            peak_index_list.append(peak_index)
        peak_index_list = np.array(peak_index_list)
        if len(peak_index_list) == 0:
            return channle_data
        # 波峰和前一个波谷的差的绝对值，以及波峰和后一个波谷的差的绝对值，放在一起排序，组成pks
        pks = np.abs(diff_channel_data[peak_index_list] - diff_channel_data[valley_index_list[0:-1]])
        pks = np.concatenate((pks, np.abs(
            diff_channel_data[peak_index_list] - diff_channel_data[valley_index_list[1:len(valley_index_list)]])))
        pks = np.sort(pks)
        if len(pks) == 0:
            return channle_data
        # 一阶差分去均值
        diff_channel_data -= np.mean(diff_channel_data)
        # 截止频率0.6，高通滤波
        diff_channel_data = filtfilt(HF_b, HF_a, diff_channel_data)
        # 通过差分信号重建ppg信号
        sig_list = [0]
        for temp_index in range(1, len(diff_channel_data)):
            sig_list.append(sig_list[-1] + diff_channel_data[temp_index])
        sig_list = np.array(sig_list)
        sig_list -= np.mean(sig_list)
        return sig_list

    @staticmethod
    def smooth_filt(data, smooth_len):
        """滑动平均滤波,去除异常值"""
        result = []
        for index in range(len(data)):
            if index < len(data) - smooth_len:
                result.append(np.mean(data[index: index + smooth_len]))
        result.insert(0, result[0])
        result.insert(0, result[0])
        result.insert(0, result[0])
        result.append(result[-1])
        result.append(result[-1])
        result.append(result[-1])
        result = np.array(result)
        return np.array(result)

    @staticmethod
    def filt_process(channle_data, FS):
        """组合滤波处理  信号滤波"""
        HF_a = [1, -3.590148110152597, 4.873443614170354, -2.963870062130241, 0.681973883426416]
        HF_b = [0.775707823645609, -3.094022833216837, 4.636644012395594, -3.094022833216837, 0.775707823645609]
        LF_a = [1, -2.174805506915435, 2.357402921272809, -1.282544135294374, 0.306064117903952]
        LF_b = [0.036936953516246, 0.031977931425171, 0.061290122632184, 0.031977931425171, 0.036936953516246]
        # 首尾拼接一秒钟采样的数据
        temp_data = np.array(list(channle_data[FS:0:-1]) + list(channle_data) + list(channle_data[-2:-2 - FS:-1]))
        temp_data = temp_data - np.mean(temp_data)
        filtered_sig = filtfilt(HF_b, HF_a, temp_data)
        filtered_sig = filtfilt(LF_b, LF_a, filtered_sig)
        filtered_sig = filtered_sig[FS:-FS]
        return filtered_sig

    @staticmethod
    def motion_filt_process(channle_data):
        """运动干扰滤波"""
        PF_b = [0.142537860520234, -0.215022543769845, 0.149319417312823, -0.215022543769845, 0.142537860520234]
        PF_a = [1, -2.67035772773509, 2.94351262911759, -1.63009314766181, 0.400438754415306]
        motion_filt_process_data = filtfilt(PF_b, PF_a, channle_data)
        return motion_filt_process_data

    @staticmethod
    def high_pass_filtering(channle_data, FS):
        """高通滤波专用方法"""
        HF_a = [1, -3.590148110152597, 4.873443614170354, -2.963870062130241, 0.681973883426416]
        HF_b = [0.775707823645609, -3.094022833216837, 4.636644012395594, -3.094022833216837, 0.775707823645609]

        # 首尾拼接一秒钟采样的数据
        temp_data = np.array(list(channle_data[FS:0:-1]) + list(channle_data) + list(channle_data[-2:-2 - FS:-1]))
        temp_data = temp_data - np.mean(temp_data)
        # 对信号进行双向滤波
        filtered_sig = filtfilt(HF_b, HF_a, temp_data)
        filtered_sig = filtered_sig[FS:-FS]
        return filtered_sig


#特征计算模块
class FeatureCalculator:
    """特征计算模块"""

    @staticmethod
    def max_peak_energy_ratio(czt_values):
        """
            最大峰能量占比
            :param czt_values:幅值列表
            :return: 返回占比值
            """
        # 获取幅值最大值
        czt_values_max = np.max(czt_values)
        #print("幅值最大值:" + str(czt_values_max))
        czt_values_sum = sum(czt_values)
        #print("幅值合:" + str(czt_values_sum))
        return czt_values_max / (sum(czt_values))

    @staticmethod
    def interfering_peaks_count(now_sig_czt, now_sig_freq):
        """
            计算干扰峰个数,要求:波峰大于0.5则为干扰峰
            :param peak_list: 波峰列表
            :return:返回干扰峰个数
            """

        # 使用 find_peaks 找到信号中的波峰索引
        value_peak_index_list, _ = find_peaks(now_sig_czt, distance=6)

        # 获取波形对应的幅值
        czt_values = now_sig_czt[value_peak_index_list]

        # 获取对应的频率的位置
        frequency_peaks = now_sig_freq[value_peak_index_list]

        # 计算peak_list中大于0.5的个数
        # 把peak_list中大于0.5的数添加进一个新列表
        peak = []
        for i in czt_values:
            if i > interference:
                peak.append(i)

        return len(peak)

    @staticmethod
    def energy_ratio(czt_values_list, czt_values):
        """
        计算算法峰能量占比
        :param czt_values_list:幅值数组
        :param czt_values: 金标所在索引位置对应的频率值 金标对应的幅值
        :return:返回算法峰能量占比
        """
        # 计算能量占比
        energy_ratio = czt_values / sum(czt_values_list)

        return energy_ratio

    @staticmethod
    def get_ac_dc_energy_ratio(window_data, Gain_list):
        """
            计算AC/DC能量比值
            :param window_data: 未处理前的数据,10s窗口的数据
            :param Gain_list: 增益列表
            :return:返回AC/DC能量比值
            """

        # 设置增益
        gain_value = np.array([10, 25, 50, 75, 100, 250, 500, 750, 1000]) / 5

        # 处理增益值
        gain_data_list = gain_value[np.array(Gain_list)]

        # 计算直流能量（DC energy）
        dc_energy = (np.array(window_data) - 2 ** 23) / gain_data_list

        # 计算交流能量（AC energy）
        x = np.array(window_data)
        # 计算信号的均值
        mean_x = np.mean(x)
        # 去掉均值
        x_prime = x - mean_x
        # 使用高通滤波
        ac_energy = SignalPreprocessor.high_pass_filtering(x_prime, 25)

        # ac/dc 能量比值列表
        ac_dc_energy_list = ac_energy / dc_energy

        # 计算 25% 和 75% 分位数
        q25 = np.percentile(ac_dc_energy_list, 25)
        q75 = np.percentile(ac_dc_energy_list, 75)
        return (q25 + q75) / 2

#峰值检测模块
class PeakDetector:
    """峰值检测模块"""

    @staticmethod
    def get_peak_of_wave(now_sig_czt, now_sig_freq):
        """
            通用波峰检测 找列表中的波峰,与波峰对应的幅值(此波峰检测修改：会影响 最大峰能量占比、算法峰检测(影响算法峰最大能力占比))   算法峰就是心率峰 或者是血氧峰,现改为统一算法峰名称
            :param now_sig_czt:
            :return: value_peak_index_list:返回波峰数组;
            now_sig_freq:频率(心率值)数组
            now_sig_czt: 波峰对应的幅值数组
            """
        # 使用 find_peaks 找到信号中的波峰索引
        value_peak_index_list, _ = find_peaks(now_sig_czt, distance=6)

        # 获取对应的频率的位置
        frequency_peaks = now_sig_freq[value_peak_index_list]
        # 获取波形对应的幅值
        czt_values = now_sig_czt[value_peak_index_list]

        new_data = list(zip(value_peak_index_list, frequency_peaks, czt_values))
        out_data = []
        # #筛选波峰,只保留大于指定幅值(谱峰)的波峰
        for peak_data, fr_data, czt_data in new_data:
            if czt_data > general_peak_threshold:  # 指定幅值(czt_)
                out_data.append((peak_data, fr_data, czt_data))
        # 将 out_data 拆分成3个列表
        peak_list, frequency_peaks_list, czt_values_list = zip(*out_data) if out_data else ([], [], [])

        return peak_list, frequency_peaks_list, czt_values_list

    @staticmethod
    def get_heart_rate_peak(now_sig_czt, now_sig_freq, now_polar_heartbeat):
        """
            金标心率匹配检测 找心率带波峰,对应的频率值,要求:最靠近心率峰与标签误差需小于10bpm,不满足这一条件时，这1s占比标记为0
            :param now_sig_czt:   z:CZT变换后归一化后的幅值 now_sig_czt
            :param now_sig_freq:   fz:CZT变换后对应的频率 now_sig_freq
            :param now_polar_heartbeat: 金标心率
            :return:返回对应最接近金标心率的波峰的 对应的幅值(czt_value),否则返回0
            """
        # 找列表中的波峰
        #value_peak_index_list, _ = find_peaks(now_sig_czt, distance=6)
        #使用通用波峰检测
        value_peak_index_list, frequency_peaks_list, czt_values_list = PeakDetector.get_peak_of_wave(now_sig_czt, now_sig_freq)
        
        # 频谱中,心率带心率的位置
        polar_heart_rate_nearest_index = \
            np.where(np.abs(now_sig_freq - now_polar_heartbeat) == np.min(np.abs(now_sig_freq - now_polar_heartbeat)))[
                0][0]
        # 寻找离心率峰最近的频谱峰的位置
        polar_hr_peaks_diff = np.abs(np.array(value_peak_index_list) - polar_heart_rate_nearest_index)
        if len(polar_hr_peaks_diff) == 0:
            nearest_polar_czt_index = 0
            return nearest_polar_czt_index
        else:
            nearest_polar_peak_index = np.where(polar_hr_peaks_diff == np.min(polar_hr_peaks_diff))[0][0]
            nearest_polar_czt_index = value_peak_index_list[nearest_polar_peak_index]
            # 获取对应的频谱峰的位置
            frequency_value = now_sig_freq[nearest_polar_czt_index]
            # 获取对应的幅值的位置
            czt_value = now_sig_czt[nearest_polar_czt_index]

            if abs(frequency_value - now_polar_heartbeat) < 10:
                return czt_value
            else:
                return 0

#文件工具类
class FileUtil:
    """文件工具类"""

    @staticmethod
    def read_numbers_from_file(file_path):
        """从文件读取数字列表"""
        numbers = []
        with open(file_path, 'r', encoding='utf-8-sig') as file:  # 使用 'utf-8-sig' 编码来自动处理BOM
            for line in file:
                line = line.strip()  # 去除换行符和空白字符
                if line:  # 确保行不为空
                    try:
                        number = int(line)
                        numbers.append(number)
                    except ValueError as e:
                        print(f"Error converting line to int: {line}. Error: {e}")
        return numbers
#频域分析模块
class FrequencyAnalyzer:
    """频域分析模块"""

    @staticmethod
    def get_czt(s, n, f1, f2, fs):
        """
        CZT变换核心方法
            s:信号
            n:CZT点数
            f1:CZT变换后最低频率
            f2:CZT变换后最高频率
            fs:采样率
            return:
            返回CZT变换后的信号和频率
            z:CZT变换后归一化后的幅值 now_sig_czt
            fz:CZT变换后对应的频率 now_sig_freq
            """
        try:
            # 确保输入信号是数值数组
            s = np.asarray(s, dtype=np.float64)
            w = np.exp(-1j * 2 * np.pi * (f2 - f1) / (n * fs))
            a = np.exp(1j * 2 * np.pi * f1 / fs)
            z = czt(s, n, w, a)
            z = np.abs(z)
            fz = (np.arange(len(z)) * (f2 - f1) / len(z)) + f1
            z = z / np.max(z)
            fz = fz * 60
            return z, fz
        except Exception as e:
            print(f"Error in get_czt: {e}")
            return None, None

    @staticmethod
    def get_data_process(window_data, fs, Gain_list=None, now_polar_heartbeat=None):
        """
                   频域特征全流程计算 (获取频域相关数据)
            :param window_data: 窗口数据(ppg光路原始数据)
            :param fs: 采样率
            :param Gain_list: 增益列表 用于计算ac/dc能量比值 默认 = None
            :param now_polar_heartbeat: 金标心率 默认 = None
            :return:
                   now_sig_czt:幅值
                   now_sig_freq:频率(心率值)数组
                   max_peak_of_wave:最大峰能量占比
                   some_peak_of_wave:干扰峰个数
                   energy_ratio_of_wave:算法峰能量占比
                   peak_list:波峰数组
                   frequency_peaks_list:波峰对应的频率数组
                   czt_values_list:波峰对应幅值数组
                   ac_dc_energy:ac/dc能量占比
            """
        # 默认算法峰能量占比返回0
        energy_ratio_of_wave = 0
        max_peak_of_wave = 0
        some_peak_of_wave = 0
        ac_dc_energy = 0
        # 1.预处理滤波等
        data = SignalPreprocessor.ppg_process(window_data, fs)
        # 2.先进行ctz转换
        now_sig_czt, now_sig_freq = FrequencyAnalyzer.get_czt(data, 256, 0.6, 4, 25)
        # 3.再进行处理
        # 找指定的大于指定幅值 波峰列表
        peak_list, frequency_peaks_list, czt_values_list = PeakDetector.get_peak_of_wave(now_sig_czt, now_sig_freq)

        # 最大峰能量占比
        max_peak_of_wave = FeatureCalculator.max_peak_energy_ratio(czt_values_list)
        # 干扰峰个数
        some_peak_of_wave = FeatureCalculator.interfering_peaks_count(now_sig_czt, now_sig_freq)
        if Gain_list is not None:
            # 获取AC/DC能量比值
            ac_dc_energy = abs(FeatureCalculator.get_ac_dc_energy_ratio(window_data, Gain_list))

        if now_polar_heartbeat is not None:
            if now_polar_heartbeat > 0:
                #获取心率波峰位置对于的幅值
                hr = PeakDetector.get_heart_rate_peak(now_sig_czt, now_sig_freq, now_polar_heartbeat)
                # print("心率波峰位置=" + str(hr))
                # 算法峰能量占比
                energy_ratio_of_wave = FeatureCalculator.energy_ratio(czt_values_list, hr)
        
        #print("最大峰能量占比=" + str(max_peak_of_wave))
        #print("干扰峰个数=" + str(some_peak_of_wave))
        #print("算法峰能量占比=" + str(energy_ratio_of_wave))
        #print("ac/dc能量比值=" + str(ac_dc_energy))
        return now_sig_czt, now_sig_freq, max_peak_of_wave, some_peak_of_wave, energy_ratio_of_wave, peak_list, frequency_peaks_list, czt_values_list, ac_dc_energy


def get_cztdata_process(window_data, fs, Gain_list=None, now_polar_heartbeat=None):
    """
               频域特征全流程计算 (获取频域相关数据)
        :param window_data: 窗口数据(ppg光路原始数据)
        :param fs: 采样率
        :param Gain_list: 增益列表 用于计算ac/dc能量比值 默认 = None
        :param now_polar_heartbeat: 金标心率 默认 = None
        :return:
               now_sig_czt:幅值
               now_sig_freq:频率(心率值)数组
               max_peak_of_wave:最大峰能量占比
               some_peak_of_wave:干扰峰个数
               energy_ratio_of_wave:算法峰能量占比
               peak_list:波峰数组
               frequency_peaks_list:波峰对应的频率数组
               czt_values_list:波峰对应幅值数组
               ac_dc_energy:ac/dc能量占比
        """
    now_sig_czt, now_sig_freq, max_peak_of_wave, some_peak_of_wave, energy_ratio_of_wave, peak_list, frequency_peaks_list, czt_values_list, ac_dc_energy= FrequencyAnalyzer.get_data_process(window_data, fs, Gain_list, now_polar_heartbeat)
    return now_sig_czt, now_sig_freq, max_peak_of_wave, some_peak_of_wave, energy_ratio_of_wave, peak_list, frequency_peaks_list, czt_values_list, ac_dc_energy